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ABSTRACT 


To  help  understand  the  stability  of  cold,  viscous  boundary  layers  in 
geophysical  contexts  such  as  lava  lakes  and  mantle  convection,  the  following 
model  problem  is  analyzed:  Beneath  a  shear- free  horizontal  boundary,  a  thin 
layer  of  very  viscous  fluid  overlies  a  deep  layer  of  less  viscous,  less  dense 
fluid.  The  initial  unstable  equilibrium  is  perturbed,  and  the  growth  of  the 
disturbance  is  followed,  including  the  nonlinear  effects  of  large  amplitude, 
by  a  long-wave  analysis.  The  result  shows  that  in  the  final  catastrophic 
growth  the  peak  thickness  of  the  upper  layer  approaches  infinity  inversely 
proportional  to  the  remaining  time.   (This  conclusion  also  applies  to  fluids 
with  power-law  rheology.)  Thus  nonlinear  effects  greatly  enhance  growth. 
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Summary  of  notation 
Pi,  ui,  P2,  M2»  AP>  g»  x»  z»  t,  u,  w  :  see  Figure  1 
P  :  local  density  ( Pi  or  P2) 
v  :  either  kinematic  viscosity 
k  :  characteristic  wavenumber 
dj,  d2  :  equilibrium  depths  of  each  layer 
6(x,y,t)  :  interface  location 
P  :  local  pressure  (Pi  or  P2) 

p  :  local  reduced  pressure  (reduced  by  pg(z-di)) 
o_   :  stress  tensor 

t^j  :  deviatoric  stress  tensor  (without  pressure  contribution) 
a  =  U2/V1  :  viscosity  ratio 
B  =  d 2/d  1  :  depth  ratio 
£  =  kdj  :  dimensionless  wavenumber 
Fij  :  reduced  force  tensor,  2D 
F(t)  :  reduced  force  for  ID  disturbance 
de(t)  :  current  "equilibrium"  thickness 
xq  :  initial  position  of  fluid  cross  section 
<5()(xo)  :  initial  thickness  profile 

a(x,t)  :  disturbance  amplitude  (dimensionless,  =  6-1) 
ao(xo)  :  initial  amplitude 

t*  :  singular  time  for  any  cross  section  to  reach  infinite  thickness 
T  -  t-t*  :  time  relative  to  singular  time 

L  :  length  of  layer,  or  wavelength  of  periodic  disturbance 
Lq  :  initial  length  or  wavelength 

n(xo,t),  f(ri),  C(n)  :  for  similarity  solutions  describing  plume 
x*(t)  :  "initial"  position  of  fluid  section  currently  at  plume  position 


a  :  arbitrary  power  for  peak  profile 

c  :  small  remainder 

°    :  linearized  growth  rate 

t(6)  =  t*-t  :  time  until  singularity  for  power-law  fluid  layer 

a(t),  A(t),  C(x,t),  f(0,  g(0  :  for  large-amplitude  similarity  solution 
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1.  INTRODUCTION 

This  work  describes  the  circumstances  under  which  the  motion  of  a 
buoyantly  unstable  horizontal  film  of  viscous  fluid  is  limited  by  normal 
viscous  stresses,  and  gives  an  analysis  showing  how  the  growth  of  disturbances 
to  the  interface  becomes  greatly  enhanced  when  the  disturbance  amplitude 
becomes  large,  leading  to  the  formation  of  plumes  in  a  finite  time.   A 
companion  work  will  show  that  in  certain  cases  this  same  dynamic  balance 
controls  the  structure  of  the  cold  boundary  layer  in  Benard  convection  with 
strongly  temperature-dependent  viscosity.   The  behavior  described  here  thus  is 
of  geophysical  relevance  to  such  flows  as  those  in  lava  lakes,  the  thermal 
convection  of  planetary  mantles,  and  possibly  also  in  convection  in  the 
Earth's  solid  core.   The  Rayleigh-Taylor  problem  considered  here  is  the 
simplest  example  of  this  dynamic  balance. 

The  problem  is  illustrated  in  figure  1.   Two  layers  of  immiscible  fluids 
are  confined  between  shear-free  boundaries,  with  the  upper  fluid  layer  being 
denser  and  much  more  viscous  than  the  lower,  and  both  fluids  are  so  viscous 
that  neither  inertia  nor  surface  tension  is  significant.   A  linearized 
analysis  (appendix  A)  shows  that  the  most  unstable  wavelength  is  long  compared 
to  the  upper  layer  thickness.   And  scaling  the  problem  shows  that  this  fastest 
growth  occurs  in  a  broad  range  of  wavelengths  for  which  the  lower  fluid  is 
effectively  passive  and  hydrostatic.   Thus  we  give  a  long-wave  analysis,  where 
the  lower  fluid  is  treated  as  inviscid,  that  exploits  the  fact  that,  for  long 
waves,  the  disturbance  amplitude  can  get  very  large  (compared  to  the  average 
layer  thickness)  while  the  slope  of  the  interface  remains  small;  this  follows 
from  mass  conservation. 
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The  results  of  this  analysis  show  the  nonlinear  effects  of  large 
amplitude.   In  particular,  it  predicts  that  thickness  maxima  in  the  layer  will 
reach  infinite  thickness  in  finite  time,  with  the  final  catastrophic  growth 
inversely  proportional  to  the  time  remaining  before  the  singular  time.   This 
conclusion  is  not  unique  to  a  Newtonian  rheology;  the  corresponding  analysis 
for  a  power-law  fluid  (appendix  B)  shows  the  same  catastrophic  inverse-time 
growth  of  large  peaks. 

The  small-slope  analysis  predicts  that  at  the  singular  time,  a  peak  will 
locally  have  the  shape  °  *  jx|'    ,  but  that  afterward,  as  the  plume  acts  as  a 
sink  of  fluid,  the  shape  will  change  to  6  a  |x|~1/2. 

Of  course  the  small-slope  approximation  must  fail  before  this,  but  we 
argue  that  the  failure  is  only  local;  the  slope  remains  small  everywhere 
except  in  an  asymptotically  small  neighborhood  around  each  developing 
singularity.   Furthermore,  these  small  regions  where  the  small-slope 
approximation  breaks  down  have  only  an  asymptotically  small  effect  en  the 
dynamics  of  the  rest  of  the  layer.   Therefore  the  small-slope  equations  will 
continue  to  apply  almost  everywhere,  even  after  the  singular  time,  when  peaks 
form  downwelling  plumes.   (There  is  a  family  of  similarity  solutions,  given  in 
appendix  C,  that  describe  plume  behavior.)   Thus  the  same  equations  describe 
the  disturbance  from  the  initial  linear  growth  through  the  rapid  large- 
amplitude  growth  to  the  final  draining  of  the  layer  by  the  plumes. 

In  the  companion  paper  on  the  analogous  thermally  driven  buoyant 
instability,  there  is  a  family  of  steady  solutions  (with  plumes),  and  we  use 
the  small-slope  equations  to  follow  the  development  from  initial  conditions 
through  to  the  final  steady  state. 
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2.  PROBLEM  STATEMENT 

Two  horizontal  layers  of  distinct  Newtonian  incompressible  fluids  are 
bounded  above  and  below  by  horizontal  shear-free  boundaries  (see  figure  1). 
The  upper  fluid  (of  density  Pi  and  viscosity  Ui)  is  denser  and  much  more 
viscous  than  the  lower  fluid  (of  density  P2  <  Pi  and  viscosity  U2  <<  Pi),  and 
the  upper  layer  is  very  thin  (di  <<  d2).   Both  fluids  are  assumed  to  be  so 
viscous  that  we  can  neglect  both  surface  tension  and  inertia.   The  latter 
assumption  requires 

ApgL3/pv2   «  i 

where  L  is  the  largest  length  scale,  Ap  =  P1-P2,  g  is  gravity,  P  is  either 
density,  and  v  is  either  kinematic  viscosity;  this  is  easily  satisfied  in 
mantle  flow. 

Initially  both  fluids  are  at  rest  in  unstable  equilibrium.   Then  at  t  =  0 
the  interface  is  slightly  disturbed;  thereafter  the  interface  position  is 
given  by  ^(x,y,t).   (Later  we  will  assume  that  the  disturbance  varies  only  in 
the  x  direction;  a  one-dimensional  disturbance  allows  some  remarkable 
simplifications,  and  illustrates  the  physical  behavior  more  clearly.) 
A  reduced  pressure  p  is  defined  by: 

P(x,y,z,t)  =  n(x,y,z,t)  +  Pg(z-di)  (2.1) 

where  p  is  the  local  density  ( p2  or  Pi).   Then  the  governing  equations  are: 

vp  =  ^v2u  (2.2a) 

v'u_  =  0  (2.2b) 

and  the  boundary  conditions  are: 
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at  z  =  0  and  at  z  =  (di+d2):  w  =  uz  =  0  (2.3a,b) 

at  z  =  6(x,y,t):  [u_]  =  0_,  (2.3c) 

[no]   =  Apg(6_dl)^5  (2.3d) 

6t  +  u<5x  +  v6y  =  w  (2.3e) 

where  LJ  indicates  the  change  in  the  enclosed  quantity  across  the  interface 
(downward),  °_   is  the  reduced  stress  tensor,  and  the  unit  normal  to  the 
interface  n  has  components: 

n  =  (  -«x,  -6y,  1  )/*/(l+6x2+6y2)  (2.4a) 

Physically,  the  conditions  (2.3a-d)  are  that  the  boundaries  are  impermeable 
and  exert  no  shear,  and  across  the  interface  both  velocity  and  stress  are 
continuous,  while  (2.3e)  is  the  kinematic  condition  that  the  interface  moves 
as  a  material  surface. 

These  equations  and  boundary  conditions,  along  with  specification  of  the 
initial  interface  position  ^(x,y,0),  define  the  problem  without  approximation. 
To  facilitate  analysis,  however,  we  limit  attention  to  the  development  while 
the  slope  of  the  interface  remains  small,  i.e.,  &x   and  ^  are  both  small. 
Then  in  the  interface  stress  conditions  we  can  approximate  the  unit  normal,  to 
lowest  order,  by: 

n  =  <  "6x>  ~V  l    }  (2'4b) 

and  neglect  any  resulting  terms  quadratic  in  &x   ancj  6  ■  it  should  be  noted 

that  the  resulting  small-slope  equations  still  retain  the  leading  nonlinear 

terms  due  to  slope  and  are  still  applied  at  z  =  <$(x,y,t),  in  contrast  to  the 

linearized  conditions  used  for  the  small-amplitude  analysis  (appendix  A). 

The  small-slope  approximation  applies  trivially  to  small-amplitude 

disturbances.   But  more  importantly,  the  small-slope  approximation  applies  to 

long-wavelength  disturbances  even  when  the  amplitude  becomes  large  (until  the 
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amplitude  is  comparable  to  the  wavelength).   Thus  we  can  use  the  simplified 
equations  to  analyze  the  large-amplitude,  nonlinear  effects  in  long-wave 
instabilities,  which  we  do  in  the  next  section. 
3.  LONG-WAVE  ANALYSIS 

The  disturbance  can  grow  in  a  variety  of  different  ways,  depending  on  the 
viscosity  ratio,  the  depth  ratio,  and  the  dimensionless  wavenumber: 

a  =    M2/P1,  B  e  d2/di,  e  e  kdi  (3.1) 

(where  the  disturbance  has  a  characteristic  wavenumber  k).   The  linearized 
analysis  in  Appendix  A  gives  the  small-amplitude  growth  rates  for  the  full 
range  of  all  parameters.   A  scaling  analysis  (see  Canright,  1987,  App.  B) 
shows  that  the  same  growth  regimes  apply  even  for  large  amplitudes,  as  long  as 
the  slope  of  the  interface  remains  small. 

Here  we  focus  on  the  case  where  the  lower  fluid  layer  is  much  less 
viscous  and  deeper  than  the  upper,  so: 

a  «  1,  3  »  1  (3.2) 

Then  there  is  a  range  of  wavelengths,  which  includes  the  fastest  growing 
wavelength,  over  which  the  growth  rate  is  nearly  constant  (see  fig  6): 

max(  a,  Aa/B)  )  «  e  «  1  (3.3) 

In  this  range,  the  growth  is  limited  by  normal  stresses  in  the  upper  fluid, 
which  moves  nearly  horizontally,  while  the  lower  fluid  is  passively  moved  by 
the  interface.   Outside  this  range,  for  waves  short  compared  to  di,  the  growth 
is  reduced  because  only  a  fraction  of  layer  1  is  mobilized,  and  for  long 
enough  waves  the  viscous  resistance  of  fluid  2  slows  the  growth. 

We  examine  the  finite-amplitude  growth  of  long-wave  disturbances  in  this 
fastest-growth  regime,  exploiting  the  fact  the  slope  of  the  interface  remains 
small  even  for  large  amplitudes  (until  the  disturbance  grows  to  the  order  of 


the  initial  wavelength).   What  follows  is  the  leading-order  asymptotic 
analysis;  since  the  growth  rate  is  constant  to  within  0(e  ),  we  cannot  expect 
this  analysis  to  predict  the  single  wavelength  giving  maximum  growth. 

For  wavelengths  in  this  range,  the  interface  motion  is  controlled  by  the 
dynamics  of  the  upper  fluid,  while  the  lower  fluid  is  effectively  inviscid  and 
hydrostatic;  the  error  we  make  in  neglecting  the  viscous  resistance  of  the 
lower  fluid  is  much  smaller  than  0(e).   Hereafter  we  refer  only  to  fluid  layer 
1  and  drop  the  subscript  1  where  clear.   The  motion  is  quasi-horizontal  as 
both  surfaces  of  the  layer  see  no  shear.   It  is  driven  by  the  horizontal 
variations  of  the  buoyant  pressure  and  resisted  by  the  normal  viscous 
stresses. 

We  now  estimate  the  scales  of  the  different  terms,  assuming  for  this 
purpose  that  the  disturbance  varies  only  in  the  x  direction.   Then  there  are 
two  length  scales,  x  ~  k" l ,  z  ~  dl,  and  5X  ~  e.   From  continuity  (and  the 
impermeable  boundary),  w/u  ~  e.   Consequently,  comparing  vertical  and 
horizontal  momentum  equations  shows  that  pz/px  ~  e,  and  the  interface  shear 
stress  condition  (2.3c)  implies  uz  ~  wx,  so  uz/ux  ~  e.   The  last  scaling  in 
turn  implies  that  the  variation  of  u  across  the  layer  is  0(e  )  smaller  than  u 


(similarly  for  p),  though  u 


zz    UXX" 


When  the  initial  disturbance  to  the  interface  varies  in  two  dimensions, 
over  a  length  scale  long  compared  to  the  thickness  of  the  upper  layer  but  not 
so  long  that  the  resistance  of  the  lower  fluid  becomes  important,  then  the 
above  scaling  still  applies.   That  is,  u,  v  and  p  are  independent  of  z,  and: 

w  =  -z(ux+Vy)  (3.4) 

all  within  an  error  of  0(e2).   From  the  normal  stress  condition,  to  0(e  ): 

"P  +  Tzz  -  APg(6-dl)  (3.5) 
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where  t^.  is  the  deviatoric  stress  tensor. 

Consider  a  force  balance  in  the  x  direction  on  a  small  column  of  layer  1 
(see  figure  10).   We  reduce  the  force  on  all  surfaces  by  the  hydrostatic 
pressure  gradient  due  to  fluid  2;  this  does  not  affect  the  balance.   Recalling 
(2.1)  and  from  (3.5),  the  difference  in  total  pressures  is: 

P1-P2  =  P  +  Apg(z-di)  =  tzz  -  Apg(6-Z)  (3.6) 

Then  the  reduced  force  balance  involves  only  the  sides  of  the  column,  since 

the  boundary  is  shear-free  and  the  reduced  force  on  the  interface  is  zero,  so: 

Ay  AX 

J       [Fxx(x+Ax)    -    Fxx(x)]    dy  +    J       [Fxy(y+Ay)    -    Fxy(y)]    dx   =   0  (3.7) 

Y  6  x 

Fxx    =    /    [-(P1-P2)    +    \x]    dz    ~~    [APg62/2   +    H\x   -    \z)\ 

°6 

F    =   /  T    67    a    f  5t   I 
xy    J      xv  u^    '   xyJ 

0 
where  Fj^  is  a  2-D  tensor  representing  the  reduced  force  in  the  layer:  any 

vertical  plane  through  the  layer  defines  a  horizontal  normal  direction,  and 

the  product  of  that  normal  and  the  tensor  F^4  gives  the  reduced  force  vector 

acting  on  that  plane.   Dividing  by  AxAy  and  taking  the  limit  as  Ax }  Ay  ■>  0: 

3/3x(Fxx)  +  3/3y(Fxy)  =  0  (3.8) 

Indeed,  a  force  balance  in  the  y  direction  shows: 

a/VjCFij)  =  0  (3.9) 

Fij  =  6("P  6ij  +  Tij) 
p  =  -Apg6/2  +  xzz 

where  P  is  the  average  (reduced)  pressure  and  °j4  is  the  Kronecker  delta. 

The  kinematic  interface  condition  becomes: 

D6/Dt  "  6t  +  u6x  +  v6y  =  -6(ux  +  vy)  (3.10) 

These  two  equations  (3.9-10)  govern  the  growth  of  long-wave  disturbances. 
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Note  that  the  above  derivation  is  independent  of  the  rheology  of  the  fluid 
layer;  any  constitutive  relation  could  be  specified. 

The  analysis  simplifies  considerably  for  a  one-dimensional  disturbance; 
we  develop  only  this  case.   In  general,  the  reduced  force  Y^a    is  a  two- 
dimensional  tensor  with  zero  divergence,  reflecting  the  balance  of  buoyant  and 
viscous  forces  in  the  layer.   However,  when  the  disturbance  is  one- 
dimensional,  then  only  one  component  of  F^  is  nonzero;  we  call  this  scalar  F. 
Then  (3.9-10)  become: 

9/8x  F  =  0  (3.11) 

D6/Dt  =  6t  +  u«x  =  -<$ux  (3.12) 

the  former  of  which  integrates  directly  to  give: 

4<5ux  +  (Apg/y)62/2  =  F(t)  (3.13) 

This  says  that  the  force  on  any  cross  section  of  the  layer,  reduced  by 
the  hydrostatic  pressure  force  due  to  the  lower  fluid,  is  the  same  everywhere 
in  the  layer,  i.e.,  F(t).   (This  result  is  not  surprising,  in  that  the  only 
outside  force  on  the  layer  is  just  the  external  hydrostatic  head.)   F(t) 
depends  on  the  conditions  specified  at  the  ends  of  the  layer,  e.g.,  if  the 
layer  had  abrupt  ends  with  no  applied  force,  surrounded  by  fluid  2,  F  would  be 
zero.   Arbitrary  end  conditions  can  thus  be  incorporated  via  F(t). 

We  will  treat  three  cases  where  F  is  simple  to  evaluate.    In  the  first 
case,  the  layer  is  of  Infinite  extent,  but  the  disturbance  is  localized,  so 
that  far  away  the  layer  remains  in  (unstable)  equilibrium,  and  F  is  constant. 
This  case  has  a  simple  analytic  solution  that  illustrates  clearly  the 
nonlinear  effects  of  finite  amplitude.   The  second  case  is  a  periodic 
disturbance  (or  equivalently ,  a  layer  with  shear-free  end  walls).   Here,  the 
requirement  that  u  =  0  at  both  ends  of  a  period  determines  F(t)  as  an  integral 
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property  of  the  shape  of  the  layer.   In  the  last  (and  simplest)  case  the  layer 
is  of  finite  extent  so  F  =  0;  this  solution  also  applies  approximately  when 
is  large  and  F  can  be  neglected. 

In  dimensionless  form,  the  governing  equations  (3.12-13)  become: 

DtyDt  =  33/3E  +  u3S/3x  =  -^u/Sx  (3.14) 

S8u73x  +  S2  =  F(t)  =  de2(t)  (3.15) 

which  can  be  combined  to  eliminate  u  by  adopting  a  Lagrangian  formulation, 
where  the  fluid  "particle"  in  this  case  is  a  material  cross  section  of  the 
layer,  to  give: 

D3/DE  =  I2   -  de2(t)  (3.16) 

Here  D/Dt  =  3/3t  +  u3/3x  is  the  time  derivative  following  a  fluid  cross 
section,  5  -    <Vdl,  t  -  at,  a  E  APgdl/(8^)  is  half  the  small-amplitude  growth 
rate,  x  =  kx,  u  -  ku/°,  k  is  the  wavenumber,  and  de(t)  =  [ F( t )/ ( APgdl 2/2) 1 l/ 2 
is  the  current  (nondimensional )  equilibrium  thickness,  i.e.,  that  thickness 
with  no  tendency  to  grow  or  shrink,  so  de(0)  =  1.   Equation  (3.16)  is  a 
Ricatti  equation,  and  can  be  integrated  analytically  for  several  different 
forms  of  de(t),  or  numerically  for  arbitrary  de(t). 

This  result  shows  that  the  vertical  motion  of  the  interface  depends  only 
on  the  local  layer  thickness  and  the  current  equilibrium  thickness  (which  may 
depend  on  the  overall  shape  of  the  whole  layer).   Thus  any  fluid  cross  section 
does  not  care  what  its  immediate  neighbors  are  doing,  and  the  growth  is 
insensitive  to  wavelength,  as  expected.   Of  course,  the  horizontal  motion  of 
any  cross  section  depends  on  the  shape  of  the  whole  layer. 
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4.  SOLUTIONS 
4.1.  LOCALIZED  DISTURBANCE 

The  most  illustrative  case  is  an  infinite  layer  with  a  disturbance  of 
finite  extent.  Then  the  force  on  the  ends  of  the  disturbed  region  remains 
constant:  de  =  1 .   This  has  a  simple  solution  (drop  tildes): 

S(xQ,t)  -  (1  +  A(xn)e2t\   A(XQ)  -  (ja.(xn)  -  1)  (4.!) 

(1  -  A(x0)e2t)  (60(x0)   +   1) 

where  6q(xq)  e  6(x0,t=0)  and  x0  identifies  a  particular  fluid  cross  section  by 

its  initial  position.   This  becomes  singular  in  finite  time:  when 

t  =  t*  =  In  |  1/ A  |  the  thickness  6  goes  to  «>  or  0,  depending  on  whether  6q  was 

greater  or  less  than  1.   (When  the  minimum  thickness  in  the  layer  first 

reaches  zero,  then  de  becomes  zero,  and  the  above  solution  becomes  invalid  for 

the  entire  layer.) 

When  the  initial  perturbation  amplitude  aQ(x0)  =  60(x0)-l  is  small,  the 

solution  reduces  to: 

6  „  L±_£n?li/2  +  O(ao2e2t)  (A. 2a) 

1  -  a0e2t/2 

and  while  the  amplitude  a(xQ,t)  e  6(x0,t)-l  remains  small: 

a  *   a0e2t  (i  +  a0(e2t_i)/2)  (4.2b) 

Initially  this  gives  the  exponential  growth  of  the  linearized  solution,  and  as 

nonlinearities  become  important  the  perturbation  growth  accelerates  for  peaks 

(a>0)  and  retards  for  troughs  (a<0).   As  a  result  of  this  and  volume 

conservation,  the  peaks  get  narrower  and  sharper  while  the  troughs  broaden  and 

flatten.   We  can  roughly  estimate  the  duration  of  the  linear  behavior  by  when 

this  small-amplitude  solution  gives  a  ~  1:  At  ~  ln|l/a0|. 
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As  the  amplitude  gets  large,  the  growth  becomes  algebraic  in  the  time 
remaining  before  the  thickness  becomes  singular  (°°  or  0): 

60  >  1:  <5  -  l/(t*-t)  (4.3a) 

60  <  1:  6  *  t*-t  (A. 3b) 

which  shows  that  the  rapid  nonlinear  growth  occurs  over  a  time  At  ~  1 ,  or, 
dimensionally,  8y/Apgd,. 

The  catastrophic  growth  shown  by  the  inverse  relation  (A. 3a)  between  peak 
thickness  and  remaining  time  is  strikingly  different  from  the  exponential 
growth  of  small  amplitudes.   For  large  peaks  the  relevant  time  scale  is  the 
time  remaining  before  the  singularity,  which  is  inversely  proportional  to  the 
current  dimensional  peak  thickness  6max:  t  ~  u/ApgSmax.   This  shows  that 
large-amplitude  effects  drastically  enhance  the  growth  of  peaks. 

To  see  the  shape  of  the  layer  requires  Eulerian  coordinates.   As  x  is  the 
integral  of  the  strain  dx/dxQ,  and  because  volume  in  the  layer  (and  in  each 

fluid  cross  section)  is  conserved,  dx/dx0  =  60/6: 

x0 
x(x0,t)  =  /   60(x0)/6(x0,t)  dx0  (A. A) 

0 

where  the  x  origin  is  chosen  at  some  stationary  fluid  cross  section. 

We  find  that,  at  least  initially,  the  overall  strain  increases,  i.e.,  the 

perturbed  section  stretches  out  in  the  x  direction,  for  any  initial 

perturbation  with  a  zero  mean.   To  see  this,  we  take  the  time  derivative  of 

the  position  L  of  the  end  of  the  disturbed  region: 

Lo 
DL/Dt  =  /   60(i-62)/62  dxQ  (4.5a) 

0 

Then  at  t  =  0,  in  terms  of  the  amplitude  a(x,t): 

L  L 

DL/Dt  =  /  (-2a  +  a2/fi)  dx  =  f   a2/6  dx  >  0  (A. 5b) 

since  /a  dx  =  0.   That  the  disturbed  region  stretches  out  in  the  x  direction 
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is  surprising,  but  is  entirely  consistent  with  the  assumption  that  the  force 
on  the  ends  remains  constant. 

For  example,  if  the  initial  perturbation  is  sinusoidal  we  get: 

6g  =  1  -  b  cos  x0  (4.6a) 

6(x0,t)  =  2  -  b  (l+e2t)  cos  Xn  U6b) 

2  -  b  (l-e2t)  cos  x0 

x(x0,t)  =  bC  sin  x0  +  (D-C)x0 

2CD       .    /(E2-b2)  tan(xD/2)  ,.    ,  , 

+  arctan  ~ - — ^ — -  (4.6c) 

/(E2-b2)  E  -  b 

C  =  (e2t-l)/(e2t+l),  D  =  4e2t/(e2t+l)2,  E  E  2/(e2t+l)       (4.6def) 
and  the  overall  strain  is  D-C  +  CD//(E2-b2) ,  which  increases  monotonically 
with  time.   This  solution  is  graphed  as  6(x)  for  various  times  in  figure  2. 
When  the  initial  amplitude  b  is  small,  the  solution  simplifies: 

6  «  1  +  A  C0S  X"  ,  A(t)  =  b  e2V2  (4.6g) 

1  -  A  cos  Xq 

x  -  - tan~l[  /(1'A)  tan(x0/2)  ]  -  x0  (4.6h) 

/(1-A2)        /(1+A) 

4.2.  CONSTANT  WAVELENGTH  DISTURBANCE 

For  the  second  case,  we  consider  the  more  physically  reasonable  situation 

of  a  layer  of  finite  extent  bounded  by  fixed  (shear-free)  end  walls  (or  the 

equivalent  situation  of  a  periodic  disturbance  of  infinite  extent).   Then  the 

end  forces  must  vary  with  time  to  give  u  =  0  at  both  ends.   Integrating  (3.9) 

in  x  shows  that  fixed  ends  require: 

L         L 
de2(t)  -  (/  6  dx)  /  (/  1/6  dx) 

0         0  K*'n 

where  L  is  the  entire  length  of  the  layer  (or  one  wavelength),  so  the 
numerator  is  just  the  total  volume  of  fluid  in  the  layer,  a  constant.   Thus 
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the  thinnest  parts  of  the  layer  have  the  greatest  effect  on  de(t).   We  find 

that,  as  a  result  of  (A. 7),  de(t)  decreases  in  general,  although  the  effect  is 

small  for  small  amplitude.   As  the  disturbance  grows  and  the  troughs  deepen 

and  broaden  while  the  peaks  become  narrower,  de  can  only  decrease  with  time. 

Combining  (4.7)  and  (3.10)  (using  dx  =  (60/6)dx0)  gives  a  single  partial 

integro-dif ferential  equation  that  can  be  solved  numerically  for  6(xQ,t)  given 

an  arbitrary  initial  profile  60(xQ): 

L  L 

D6/Dt  =  62  "  (  f  60  dx0)  /  (f  60/<S2  dxQ)  (4.8) 

For  large  amplitude,  the  main  effect  of  fixed  wavelength  is  to  slow  the 
growth  of  troughs;  in  fact,  6  is  prevented  from  reaching  zero.   In  contrast, 
peak  growth  is  only  slightly  accelerated.   When  peaks  become  large,  then  de 
becomes  negligible  in  comparison,  as  in  the  constant  de  case,  but  sooner  here 
since  de  decreases.   Thus  large  peaks  show  the  same  catastrophic  growth  to 
infinite  thickness,  given  by  (4.3a).   Therefore  the  growth  of  large  peaks, 
where  de  is  unimportant,  is  insensitive  to  (reasonable)  end  conditions,  and 
our  previous  conclusion,  that  the  growth  of  peaks  is  dramatically  enhanced  by 
large-amplitude  effects,  applies  in  general. 

Qualitatively  this  case  is  very  similar  to  the  previous,  as  shown  by  the 
profiles  in  figure  3,  except  that  the  wavelength  remains  constant.   Figure  4 
compares  the  growth  of  the  peak  for  fixed  wavelength  with  that  for  constant 
end  forces,  each  for  an  initial  sinusoidal  disturbance. 
4.3.  LARGE-AMPLITUDE  BEHAVIOR 

Where  6  >>  de ,  we  can  approximate  (3.16)  by 

D6/Dt  =  62  (4.9) 

(This  would  apply  without  approximation  if  the  layer  were  of  finite  extent, 
surrounded  by  fluid  2.)   This  integrates  to: 
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6(x0,t)  =  l/(l/60(x0)  -  t)  (4.10) 

where  again  6q(xq)  is  the  initial  profile.   Thus  each  cross  section  reaches 
infinite  thickness  at  time  t*  =  1/6q;  this  is  the  same  growth  as  (4.3a). 
For  this  case,  (3.15)  simplifies  to: 

6  +  ux  =  0  (4.11) 

so 

X 

u(x,t)  =  -  /  6(s,t)  ds  (4.12) 

0 
where  the  origin  is  assumed  to  be  stationary.   As  long  as  6  remains  finite 

everywhere,  we  can  use  6  dx  =  60  dxQ,  so: 

x0 
u(x0)  =  -  /   60(s)  ds  (4.13) 

0 

which  shows  that  each  cross  section  moves  at  a  uniform  speed  until  the  maximum 

thickness  reaches  the  singular  time.   These  results  (4.10,4.12)  apply  locally 

around  any  large  peak  where  de  is  negligible. 

5.  PLUME  FORMATION 

Here  we  examine  what  happens  to  a  thickness  maximum  as  it  approaches  the 
singular  time  (t*).   We  will  show  how  the  peak  becomes  a  plume,  where  the 
dense,  viscous  fluid  drains  down.   We  argue  that  the  small-slope  equations 
still  describe  the  behavior  of  the  fluid,  except  in  an  asymptotically  small 
neighborhood  of  the  plume.   This  is  possible  because  that  neighborhood  has 
only  an  asymptotically  small  effect  on  the  stress  in  the  layer.   Thus  we  can 
describe  the  shape  of  the  layer  around  the  plume,  and  the  small-slope 
equations  can  be  used  all  the  way  from  initial  conditions  through  the  final 
draining  of  the  layer. 

As  the  maximum  grows,  clearly  the  small-slope  assumption  must  break  down 
locally  before  the  peak  reaches  infinite  thickness;  at  what  thickness  it 
breaks  down  depends  on  the  initial  wavelength  of  the  disturbance.   Since  the 
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initial  disturbance  wavelength  is  asymptotically  long,  or  equivalently ,  the 
layer  is  initially  asymptotically  thin  (e  =   kdj  <<  1),  then  when  the  peak 
reaches  infinite  thickness,  that  portion  of  the  layer  around  the  peak  where 
the  physical  slope  of  the  interface  is  0(1)  or  greater  will  be  asymptotically 
narrow  compared  to  the  whole  wavelength.   This  follows  from  mass  conservation, 
and  from  the  shape  near  the  peak  approaching  (as  we  will  show)  an  integrable 
negative  power  of  x  (where  for  convenience  we  have  chosen  x  =  0  at  the  peak). 
This  point  is  illustrated  in  figure  5.   So  the  small-slope  approximation 
continues  to  apply  almost  everywhere  in  the  layer;  what  is  needed  is  a 
description  of  the  effect  of  the  peak  or  plume  on  the  rest  of  the  layer 

In  a  companion  paper,  we  give  a  large-slope  analysis  appropriate  to  the 
region  around  a  peak  where  the  small-slope  approximation  no  longer  applies. 
There  the  flow  is  extensional,  nearly  vertical,  driven  by  buoyancy  and  limited 
by  normal  viscous  stresses.   We  find  that  large-slope  effects  do  not  slow  down 
the  growth,  they  only  affect  the  detailed  shape  of  the  peak.   Specifically,  we 
find  that  the  catastrophic  behavior  described  by  (4.3a)  still  applies  (except 
for  a  numerical  coefficient  close  to  1).   Physically,  there  is  nothing  to 
prevent  the  fluid  from  flowing  down.   In  this  way,  the  peak  extends  to  become 
a  plume,  where  the  viscous  fluid  continues  to  drain  down.   Of  course,  at  some 
point  the  fluid  will  reach  the  lower  boundary  or  the  extending  peak  will 
become  so  long  that  the  viscous  resistance  of  fluid  2  becomes  important.   In 
the  former  case,  only  the  details  of  the  plume  shape  are  changed,  but  in  the 
latter,  the  driven  flow  In  fluid  2  could  affect  the  dynamics  of  the  upper 
layer. 

To  the  rest  of  the  layer  the  plume  appears  as  an  isolated  singularity,  a 
sink  of  fluid.   The  horizontal  force  balance  must  still  apply  even  to  the 
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plume,  and  so  the  reduced  force  in  the  layer  F(t)  is  continuous  across  the 
plume.   Recall  that  F(t)  is  an  integral  property  of  the  whole  layer;  for  a 
fixed-wavelength  disturbance  the  dimensionless  form  is  given  by  (4.7). 
Because  the  region  where  the  small-slope  approximation  fails  is  integrable  (as 
it  must  be,  since  no  new  mass  is  created)  and  asymptotically  thin,  its  effect 
on  F(t)  is  negligible.   Equation  (4.7)  shows  that  where  6  is  largest  has  the 
smallest  effect  on  F.   (Even  after  fluid  starts  draining  out  of  the  layer,  the 
numerator  in  (4.7)  still  just  gives  the  volume  of  the  layer,  regardless  of  its 
shape. ) 

As  an  example,  consider  what  happens  near  the  maximum  when  the  initial 
disturbance  is  a  small-amplitude  sinusoid.   (For  simplicity,  we  assume  the 
constant-force  end  conditions,  but  as  F  has  little  effect  on  a  large  peak  the 
results  apply  to  more  general  end  conditions.)   The  previous  solution  (4.6g,h) 
shows  that  the  peak  becomes  singular  as  A  =   be2t/2  approaches  unity.   Then  in 
(4.6h)  the  argument  of  the  inverse  tangent  approaches  zero  (for  xQ  away  from 
tt),  so  near  the  singular  time  (4.6h)  becomes: 

x  «  2tan(x0/2)  -  x0  (5.1a) 

for  (1-A)  «  1,  and  (1-A)1/3  <<  (tt-xq).   Then  near  the  peak  (x0  «  1): 

x  m   xQ3/i2  (5.1b) 

Setting  A  =  1  in  (4.6g)  gives  the  thickness  profile  at  the  singular  time: 

5  -  cot2(x0/2)  (5.2a) 
and  near  the  peak: 

6  «  (x0/2)~2  m   (3/2  x)"2/3  (5.2b) 
This  shows  that  at  the  singular  time  the  peak  becomes  proportional  to  an 
integrable  negative  power  of  x. 
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In  fact,  any  smooth  initial  peak  results  in  the  same  power  of  x.   For 
x0  <<  1,  the  smooth  (symmetric)  initial  peak  is  asymptotically: 

60  •  1  +  b(l  -  cxQ2/2)  +  OCbxQ4)  (5.3) 

where  b  is  the  (small)  initial  amplitude  and  c  is  some  positive  constant. 
Then  from  (4.2a): 

.    1  +  A(l  -  c  xn2/2) 

6  *  u (5.4a) 

1  -  A(l  -  c  xQ2/2) 
where  again  A(t)  =  be2t/2.   Near  the  singular  time,  A  -►  1  and: 

6  + (5.4b) 


(1-A)  +  c  x02/2 
From  continuity,  dx/dxQ  =  60/6  «  1/6  (for  small  initial  amplitude),  so  for 
Xq  <<  1  and  as  A  •>  1 : 

x  +  (l-A)x0/2  +  c  x03/12  (5.5) 

At  the  singular  time,  A  =  1,  and  near  the  singularity: 

6  «  4/(c  x02)  M  c-1/3  (3x/2)"2/3  (5.6) 

For  the  sinusoidal  disturbance,  c  =  1  and  (5.6)  reduces  to  (5.2b). 

As  the  singular  time  is  approached,  the  peak  can  be  described  by  a 
similarity  solution.   Defining  n(*0,t)  by: 

(1-A(t))n2  =   c  x02/2  (5.7a) 

reduces  (5. 4b, 5. 5)  to: 

f(n)  =    (l-A)6/2  =  l/(l+n2)  (5.7b) 

£(n)  =  3/(2c)(l-A)"3/2  x  =  n3  +  3n  (5.7c) 

This  is  a  particular  case  of  the  general  similarity  solution  for  large  6  given 
in  appendix  D.   The  comparison  is  clarified  by  noting  that: 

A  e  e2i,  T  =  t-t*  (5.8a) 

so  that  as  A  ■>  1  ,  t  -»■  0 ,  then  to  leading  order  in  t: 
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(1-A)  +  2t  (5.8b) 

The  general  similarity  solution  shows  that  a  peak  of  the  more  general  form 
60  "  1  +  b(l  -  c | Xq I n) ,  where  n  >  0,  gives  a  singularity  of  the  form 
S(x,t*)  <r  |x|-m,  where  m  =  n/(n+l).   In  other  words,  any  (reasonable)  initial 
peak  becomes  ,  at  t*,  a  singularity  proportional  to  an  integrable  negative 
power  of  x. 

After  the  peak  becomes  singular  (i.e.,  a  plume)  the  shape  of  the  layer 
around  the  plume  can  be  calculated  using  the  large-amplitude  forms  of  the 
equations  given  in  section  4.3.   For  simplicity,  we  assume  6  is  symmetric  in 
x,  decreasing  monotonically  away  from  the  origin,  and  we  only  consider  x  >   0. 
Also,  we  choose  the  time  origin  such  that  the  "initial"  amplitude  60  >>  de  in 
the  region  of  interest.   Then  the  large-amplitude  equations  apply: 

6(x0,t)  =  l/(l/60(x0)  -  t)  (4.10) 

x 
u(x,t)  =  -  /  6(s,t)  ds  (4.12) 

0 
The  peak  becomes  infinite  at  t*  =  1/ 60(0) .   Thereafter,  the  layer  must 

move  in  such  a  way  that  each  fluid  cross  section  reaches  the  origin  at  the 

same  time  that  it  reaches  infinite  thickness.   In  other  words,  if  for  t  >   t* 

we  designate  the  fluid  element  xQ  at  the  singularity  by  xQ  =  x*(t),  so 

x*(t*)  =  0,  then  6o(x*(t))  =  l/t>  and  since  <$0  is  an  invertable  function  by 

the  assumptions  above,  this  uniquely  defines  x*,  and  thus  also  determines  the 

strength  ( 60(x*)dx*/dt )  of  the  sink  at  the  singularity,  as  a  function  of  time. 

With  these  large-amplitude  equations  we  can  determine  the  shape  around  the 

singularity  for  all  time. 

Now  we  consider  the  fate  of  the  inverse  power  of  x  singularity  that  forms 

at  t*: 

at  t  =  0:  60(x0)  =  5(x=x0,t=t*)  =  b  x0-ct,  0  <  a  <  1  (5.9a) 
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where  we  have  chosen  the  reference  time  at  t*,  when  T  =  t  -  t*  =  0,  and  b  is 
some  positive  amplitude.   Then  the  fluid  cross  section  xc  at  the  origin  is 
specified  by  x0  =  x*(x),  where: 

x*(T)  =  (bi)1^  (5.9b) 

Using  (A. 10)  shows  that: 

6(x0,T)  =  b/(xQa  -  x*a)  (5.9c) 

so: 

x0 
x(x0,T)  =  /    60(s)/6(s,T)  ds 

x* 

=  (x0  -  x*)  -  (x*a/(l-a))(x0l-a  -  x*l_a)  (5.9d) 

Far  away  from  the  plume  (xQ  >>  x*),  the  profile  is  still  the  starting  profile 
from  t  =  0  (x  «  Xq,  6  *  bx~c).   However,  very  close  to  the  plume: 

x0  =  x*(l+£),  e   «    1  (5.9e) 

x  «  ox-Ae2/2  (5.9f) 

6/b  *  l/(aex*a)  «  (bx) l  ><  2a)~1//(  2ax)  (5.9g) 

This  shows  that  asymptotically  near  the  plume,  the  singularity  goes  like 
l//x,  with  a  scale  that  varies  in  time.   This  asymptotic  x  dependence  is 
actually  independent  of  the  starting  conditions  (as  shown  by  Canright,  1987, 
App.  C).   From  (5.9g)  it  is  clear  that  whether  6  near  the  plume  grows  or 
shrinks  is  determined  by  whether  a  <  1/2  or  a  >  1/2,  respectively.   The 
special  case  a  =  1/2  gives  a  steady  solution.   For  a  >  1/2,  the  fluid  drains 
away  down  the  plume  faster  than  it  comes  in  from  the  sides,  and  the  square- 
root  singularity  diminishes  with  time  as  it  spreads  out,  to  match  onto  the 
nearly  undisturbed  profile  x~a.   This  would  be  the  eventual  fate  of  an 
initially  (t  =  0)  smooth  maximum,  which  gives  a  =  2/3.   Conversely,  if 
a  <  1/2,  the  square-root  singularity  grows  as  it  spreads,  fed  from  the  sides 
faster  than  it  can  drain  fluid  away.   (To  get  a  <  1/2  would  require  a  cusp- 
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like  initial  maximum,  which  may  not  be  physically  realistic.)   This  solution 
(5.9a-g)  is  again  a  particular  case  of  the  general  large-amplitude  similarity 
solution  of  appendix  D.   To  see  this,  define  ti(x0,t): 

n(x0,i)  e  x0/x*(T),  n  >    1  (5.10a) 

f(n)  =  t6  =  l/(na  -  1)  (5.10b) 

5<n)  E  x/x*(T)  =  n  -  (n1_a-a)/(l-a)  (5.10c) 

With  the  above  description  of  how  a  plume  first  forms  and  how  it  behaves 
afterward,  the  small-slope  equations  can  be  used  to  follow  the  development  of 
the  instability  from  initial  conditions  through  rapid  large-amplitude  growth 
all  the  way  to  the  draining  away  of  the  fluid  down  the  plumes.   The  results 
will  be  inaccurate  wherever  the  physical  slope  of  the  interface  is  not  small, 
but  such  regions  comprise  only  a  small  fraction  of  the  domain.   The  only 
assumption  is  that  a  plume  does  not  exert  any  net  horizontal  force  on  the 
surrounding  layer.   This  assumption  may  break  down  if  a  plume's  length  becomes 
so  much  greater  than  the  initial  wavelength  that  the  flow  it  drives  in  fluid  2 
becomes  dynamically  significant. 
6.  CONCLUSIONS 

The  central  concern  of  this  work  is  the  nonlinear  interactions  between 
buoyant  forces  and  normal  viscous  stresses  that  occur  in  a  buoyantly  unstable 
viscous  layer  under  a  shear-free  horizontal  boundary  and  over  a  much  less 
viscous  fluid,  for  long-wave  disturbances  in  the  range  where  the  lower  fluid 
is  dynamically  unimportant.   After  the  initially  uniform  thickness  of  the 
layer  has  been  perturbed  slightly,  the  early  growth  of  the  (small) 
perturbation  is  exponential;  the  perturbation  keeps  its  shape  while  it  grows. 
But  when  the  nonlinearities  become  important,  the  thicker  parts  of  the  layer 
thicken  more  rapidly  while  the  thinner  parts  thin  more  slowly,  giving  in 
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general  sharp  peaks  with  broad,  flat  troughs  in  between.   The  accelerating 
growth  of  peaks  leads  to  infinite  thickness  at  some  time  t*  (which  depends  on 
initial  conditions),  and  the  final  catastrophic  growth  of  the  peak  thickness  6 
is  algebraic:  5  =  u/Apg(t*-t).   (In  fact,  this  same  sort  of  inverse-time 
catastrophic  growth  is  also  predicted  for  a  power-law  fluid,  though  the 
coefficient  is  different.)   This  shows  how  large-amplitude  growth  is 
fundamentally  different  from  small-amplitude  growth;  large-amplitude  effects 
dramatically  enhance  the  growth  of  peaks. 

These  are  the  results  of  a  small-slope  analysis;  of  course,  where  the 
peaks  become  very  large  the  small-slope  approximation  breaks  down  and  there 
the  flow  becomes  fully  two-dimensional.   But  as  a  companion  work  that  employs 
a  large-slope  analysis  will  show,  the  growth  of  the  disturbance  at  large 
slopes  is  essentially  the  same  as  that  predicted  here,  with  catastrophic 
inverse-time  peak  growth;  although  the  details  of  the  peak  shape  do  change, 
this  only  changes  the  prediction  by  a  numerical  factor  of  order  1. 

Furthermore,  the  small-slope  equations  will  continue  to  apply  to  the 
layer  even  after  the  formation  of  downwelling  plumes,  except  in  an 
asymptotically  narrow  neighborhood  around  each  plume.   This  is  possible 
because  the  plumes  do  not  change  the  horizontal  force  balance  (unless  the  flow 
they  drive  in  the  lower  fluid  becomes  dynamically  significant).   Applying  the 
equations  up  to  the  singular  time  shows  that  the  first  singularity  should  have 
the  local  shape  6  <*  |x|~2/3,  but  that  afterwards,  as  the  plume  drains  the 
layer,  the  singularity  changes  shape  to  6  «  |x|~1/2.   This  behavior  is 
clarified  by  a  family  of  similarity  solutions,  appropriate  where  6  is  large. 

Thus  the  small-slope  analysis  can  be  extended  to  describe  the  development 
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of  the  layer  for  all  time,  except  in  those  narrow  regions  where  the  physical 
slope  of  the  interface  becomes  large. 
APPENDICES 
A.  LINEARIZED  SOLUTION 

Here  we  calculate  initial  growth  rates  for  the  general  problem  with 
arbitrary  wavelength,  depths,  and  viscosities.   Initially  the  perturbations 
are  small,  so  the  problem  can  be  linearized  by  applying  the  interface 
conditions  at  z  =  dj  and  assuming  the  slope  remains  small,  but  of  course 
including  the  buoyant  pressure  term  due  to  the  perturbations.   (This  is  an 
extension  of  the  analysis  given  by  Whitehead  and  Luther,  1975,  to  include  the 
effects  of  finite  depth  in  fluid  2.)   Thus  the  interface  conditions  (2.3c-f) 
become : 

at  z  =  dj:  [u]  =  [u(uz+wx)]  =  0  (A. la) 

[-p  +  2uwz]  =  Apg(6-d1)  (A. lb) 

6t  =  w  (A.lc) 

where  again    indicates  the  jump  in  value  from  fluid  1  to  fluid  2.   Then  the 
equations  are  separable,  and  a  simple  analytical  solution  is  possible.   The 
perturbation  is  assumed  to  be  sinusoidal,  but  as  the  problem  is  linear,  any 
(small)  initial  condition  can  be  constructed  by  superposition. 

The  solution  is  given  below,  where  the  subscripts  '2  and  1  refer  to  the 
two  fluids,  k  is  the  wave  number,  a(t)  is  the  dimensionless  amplitude, 
Z  e  z  -  (dj+d2)  is  the  coordinate  in  fluid  2,  a  =  U2/yi  *s  tne  viscosity 
ratio,  and  in  the  coefficients  A,  B,  E,  F,  and  their  common  denominator  D 
these  abbreviations  are  used  :  k   e  2kd1?  K  =  2kd2,  c  e  cosh(k),  C  e  cosh(K), 
s  e  sinh(k),  S  e  sinh(K). 


23 
6   -  dj  [1  +  a(t)  cos(kx)]  (A.2a-j) 

Vj  =  (Apgd!/Ulk2)  a(t)/2  sin(kx)  [A  sinh(kz)  +  B  kz  cosh(kz)] 
Y2  -  (Apgdi/pjk2)  a(t)/2  sin(kx)  [E  sinh(kZ)  +  F  kZ  cosh(kZ)] 
Pi  =  (Apgdj)  {  -a(t)  B  cos(kx)  cosh(kz)  } 
P2  "  (Apgdx)  {  -a(t)  aF  cos(kx)  cosh(kZ)  } 
A  E  (-1/D)  {  [2(S-K)  +  ak(C-l)]  sinh(kd1)  + 

[k(S-K)  +  a(kK  +  2(C-1))]  coshd^)  } 
B  e  (2/D)  {  (S-K  +  aK)  sinhCkd^  +  a(C-l)  coshdcdj)  } 
E  =  (1/D)  {  [2a(s-k)  +  K(c-l)]  sinh(kd2)  + 

[aK(s-k)  +  Kk  +  2(c-l)]  cosh(kd2)  } 
F  e  (-2/D)  {  [a(s-k)  +  k]  sinh(kd2)  +  (c-1)  cosh(kd2)  } 
D  e  (S-K)(s+k)  +  2a(Cc-l+Kk)  +  a2(S+K)(s-k) 
Then  from  6t  =  w(z=d1)  we  get  the  growth  rate: 

a(t)  =  a(0)  eat  (A. 3a) 

a  =   (Apgdj/yj)  o  (A. 3b) 

1   (S-K)(c-1)  +  q(s-k)(C-l) (.    .  v 

a  =  —  \h.5c) 

k   (S-K)(s+k)  +  2a(Cc-l+Kk)  +  a2(S+K)(s-ic) 

The  symmetry  of  the  problem  is  apparent  in  the  solution.   Had  we  chosen 
fluid  2  as  the  reference  rather  than  fluid  1,  the  solution  would  have  the  same 
form.   Thus  we  can,  without  loss  of  generality,  assume  that  fluid  2  is  the 
deeper  layer:  K  >   k. 

This  solution  is  governed  by  three  independent  dimensionless  parameters: 
the  non-dimensional  wavenumbers,  or  depths,  k  and  K  (or  equivalently  k  and  the 
depth  ratio  6  =   d2/di),  and  the  viscosity  ratio  a  =  y2/u1.   It  should  be  noted 
that  a  is  a  monotonically  decreasing  function  of  a,  i.e.,  if  we  increase  y2 
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while  iij  stays  constant  the  growth  rate  can  only  decrease.   Also  note  that  if 
6  ♦  »  then  (A. 3c)  reduces  to: 

5-1   (c-»  +  °(S-g)  .  (A. 3d) 

ic   (s+ic)  +  2ac  +  a2(s-k) 

which  is  just  the  result  of  Whitehead  and  Luther  (1975).   Figure  6a  shows  the 

effects  of  viscosity  ratio  on  c(k)  for  8  =  °°;  figure  6b  shows  the  effects  of 

finite  depth  for  a  «   1 . 

When  6  >>  1  there  are  well  defined  regimes  of  growth  where  different 
force  balances  are  dominant.   These  are  shown  schematically,  with  the 
corresponding  growth  rates,  in  figure  7.   As  we  discuss  the  various  growth 
mechanisms  below,  it  should  be  kept  in  mind  that  the  same  mechanisms  continue 
to  apply  as  long  as  the  slope  of  the  interface  remains  small,  which  for  long 
waves  includes  large-amplitude  growth.   (This  point  is  supported  by  the 
scaling  argument  in  Canright,  1987,  App.  B.) 

For  sufficiently  short  waves  (k  >>  min( 1 ,max( B"1 , a"1 / 3) ) )  the  disturbance 
sees  neither  boundary  and  ~a   +   [k(l+a)]_1,  so  if  one  viscosity  is  much  larger, 
that  one  limits  the  growth.   The  dominant  mechanism  here  is  that  the  vertical 
motion  of  the  interface  is  resisted  by  normal  viscous  stresses  in  the  more 
viscous  fluid,  and  the  growth  rate  diminishes  with  decreasing  wavelength. 

At  the  other  extreme,  for  sufficiently  long  waves 
(k  <<  min(6-1 ,/( 6/a) ,/(a8) ))  the  boundaries  confine  the  flow  to  be  mainly 
horizontal,  limited  by  shear  stresses  at  the  interface.   Then 
a  -►  k2(l+g/a)/l2,  and  the  controlling  viscosity  depends  on  if  a  >  g  or  not. 

Between  these  extremes,  the  waves  are  long  compared  to  layer  1  so  the 
motion  of  the  interface  is  primarily  horizontal,  and  there  are  four  different 
regimes.   In  one  (a-1  «  k  «  a~1/3,  1  «  a  <<  B3),  fluid  2  is  relatively 
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immobile  and  the  growth  is  limited  by  the  shear  across  layer  1,  giving  the 
same  growth  rate  as  the  previous  case,  i.e.,  k2/12.   (This  possibility  was 
apparently  overlooked  by  Whitehead  and  Luther.)   In  another 
(B-1  <<  k  <<  min(a,a-1)),  the  slight  resistance  of  fluid  2  (the  rate- 
controlling  viscosity  is  y2)  gives  a  small  shear  gradient  across  layer  1, 
which  over  the  long  wavelength  is  sufficient  to  balance  the  buoyancy,  giving 
a  •*■  k/4a. 

In  the  other  two  regimes,  the  less  viscous  fluid  is  effectively  passive, 
and  the  primarily  horizontal  motion  of  the  more  viscous  layer  is  limited  by 
normal  viscous  stresses.   The  resulting  growth  rate  is  nearly  independent  of 
wavelength,  and  includes  the  maximum  growth  rate  possible  for  a  given 
viscosity  contrast  a  where  this  behavior  occurs.   (For  other  a,  i.e., 
1  <<  a  <<  B3»  we  expect  the  fastest  growth  at  the  crossover  between  short-  and 
long-wave  behavior,  at  kmax  *  2.9/a1/3,  Omax  a   0.23/a2/3.)   When  fluid  2  is 
much  more  viscous  (a  >>  B3)>  this  regime  (Ab/cx)  <<  k  <<  B~ 1 )  gives  a  -»■  BMa, 
while  for  fluid  1  more  viscous  ( a   <<  D  this  regime  (max(a, A aB) )  <<  k  <<  1) 
gives  a  +   1/4. 

In  the  last  case,  consideration  of  higher-order  effects  shows  that,  for 
B"5  «  o  «  1,  a  -  ( 1/4)  ( 1  -  (a/k  +  kV720)).   This  broad  maximum  peaks  at 
^max  *  (180a)1/5  =  2.8  a1/5  and  omax  «  ( 1/4)  ( 1  -  0.44  ctk/b).      As  an  indication 
of  the  flatness  of  this  peak,  for  the  range  ak/5   <  k  <  5.2  a1/2°,  a  >  ( 1/4)  ( 1 
-  a1/5).   When  a  <<  B-5,  the  finite  depth  modifies  the  maximum  growth  rate 
giving  o  «  ( 1/4) ( 1  -  (3a/Bk2  +  k^/720))  with  a  broad  peak  at 
kmax  *  3.2(a/B)1/6  and  omax  -  (1/4)(1  -  0.44(a/B)2/ 3 ) . 

The  present  work  is  only  concerned  with  the  case  of  a  thin,  viscous  layer 
over  a  less  viscous,  deep  layer,  so  a  <<  1  and  B  >>  1 .   We  further  restrict 
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our  consideration  to  the  mechanism  giving  the  fastest  growth,  i.e.,  the  last 
regime  considered  above,  where  a   ~   1/4.   The  lowest-order  finite-amplitude 
analysis  (section  3)  therefore  predicts  that  the  growth  is  independent  of 
wavelength. 
B.  POWER-LAW  FLUID 

The  long-wave  quasi-one-dimensional  analysis  of  section  3  is  not  limited 
to  Newtonian  rheology;  any  constitutive  relation  can  be  accomodated.   The 
important  point  is  that  both  surfaces  of  the  layer  are  shear-free  (in  the 
wavelength  range  where  the  lower  fluid  is  passive),  so  throughout  the  layer 
shear  stresses  are  0(e)  smaller  than  normal  stresses,  and  the  latter  are 
independent  of  z  to  0(e2).   For  a  one-dimensional  disturbance  (3.9-10)  become: 
Apg62/2  +  26txx  =  F(t)  (B.l) 

D6/Dt  =  5t  +  u6x  =  "6ux  (B.2) 

where  x^j  is  the  deviatoric  stress  tensor. 

Below,  we  consider  the  particular  case  of  a  fluid  with  a  power-law 
rheology.   In  mantle  convection,  the  oceanic  lithosphere  makes  up  most  of  the 
cold,  stiff  upper  thermal  boundary  layer.   The  lithosphere  behaves  like  a 
rigid  plate  rather  than  a  viscous  fluid  layer,  so  the  results  for  a  Newtonian 
fluid  may  be  a  poor  model  for  destabilization  of  the  lithosphere  and  the 
initiation  of  subduction.   A  better  model  might  be  a  fluid  where  the  stress 
depends  on  some  power  of  the  strain  rate.   This  will  not  accurately  describe 
fracturing,  but  can  at  least  incorporate  weakening  at  high  rates  of  strain. 

For  a  power-law  fluid  In  quasi-one-dimensional  flow: 

Txx  =  yrl"x|n  (B.3) 

where  yr  has  the  appropriate  units.   Combining  (B.l-3)  we  get  a  single 
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(dimensionless)  Lagrangian  equation  describing  the  evolution  of  the  thickness 
of  the  layer  as  we  follow  a  material  cross  section: 

D6/Dt  =  sgn(6-de(t))  6  [ | 62-de 2( t ) | / 6]m  (B.4) 

where  m  =  1/n,  §  E  <s/h,  t  5  (Apgh/8Pr)m  t,  and  de(t)  is  the  nondimensionalized 
current  equilibrium  thickness  as  before.   For  the  special  case  of  a  Newtonian 
fluid  (n  =  m  =  1)  this  reduces  to  (3.10).   When  de(t)  is  specified,  (B.4)  can 
be  integrated  numerically  from  the  initial  thickness  60(xn).    (Now  drop 
tildes.) 

For  an  infinitely  long  layer  with  a  localized  disturbance  (i.e.,  constant 
end  forces:  de  =  1)  and  m  an  odd  integer: 

D6/Dt  =  (62-l)m/6m~1  (B.5) 

This  can  be  integrated  by  parts.   For  example,  if  m  =  3: 

t(6)  =  l/4[  6/(62-l)2  +  l/2[  6/(62-l)  +  1/2  In | ( 6-1 )/(6+l ) |  ]]  (B.6) 
where  t(6)  =  t*-t  is  the  time  remaining  before  the  thickness  of  the  fluid 
cross  section  goes  to  00  or  0.   Profiles  (6(x)  for  various  t)  are  shown  in 
figure  8  for  m  =  3  and  m  =  9,  from  an  initial  sinusoidal  disturbance,  using 
solutions  of  the  above  form. 

The  effect  of  increasing  m  is  seen  to  be  concentration  of  the  deformation 
into  narrow  regions  at  the  centers  of  peaks  and  troughs,  while  elsewhere  the 
profile  becomes  linear  in  x.   The  growth  remains  slow  initially,  then  suddenly 
becomes  catastrophic  after  a  certain  threshold  has  been  reached.   A  trough 
reaches  this  threshold  and  necks  off  much  sooner  than  a  peak  of  equal  initial 
amplitude  blows  up.   This  can  be  seen  in  figure  9,  which  shows  the  growth  of 
positive  and  negative  perturbations  over  time. 

When  the  amplitude  a  =  6-1  is  small,  the  approximate  solution  is: 

a  -  a0/[l  -  (m-l)(2a0)m-lt]  l/(m-l)  (B.7) 
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where  aQ  is  the  initial  amplitude.   This  slow  growth  lasts  for  a  period  of 

roughly  At  ~  1/ f (m-1 ) ( 2a0 )m_  1  ] .   For  large  amplitudes,  the  growth  again 

becomes  algebraic  in  the  time  remaining  before  the  singularity  is  reached  at 

t  =  t*: 

60  >  1:   6  -  l/(m(t*-t))  (B.8a) 

60  <  1:   6  «  m(t*-t)  (B.8b) 

This  catastrophic  growth  occurs  over  a  time  scale  At  ~  1/m  (dimensionally 

(8yr/Apgh)m/m). 

As  for  the  Newtonian  fluid,  a  disturbed  region  under  constant  end  forces 

will  tend  to  stretch  out  to  some  extent  in  the  x  direction.   To  keep  the 

wavelength  constant,  the  end  forces  must  vary  in  time  to  give: 
L 
f  [(s2-de2(t))/6]m  6  dx  =  0  (B.9) 

While  we  have  not  done  a  numerical  calculation  for  this  case,  we  speculate 
that  enforcing  constant  wavelength  does  not  significantly  alter  the 
qualitative  behavior,  except  to  inhibit  the  layer  from  necking  off.   (When  a 
fluid  cross  section  reaches  zero  thickness,  locally  the  strain  must  become 
infinite;  for  smooth  initial  conditions,  this  is  incompatible  with  constant 
overall  strain.)   Also,  peak  growth  may  be  slightly  enhanced. 

The  main  point  is  that  when  peaks  get  large  (compared  to  the  current 
equilibrium  thickness),  the  large-amplitude  effects  still  produce  catastrophic 
growth  of  the  form  6  <*  l/(t*-t),  as  for  a  Newtonian  fluid,  where  t*-t  is  the 
time  remaining  before  the  singularity. 
C.  LARGE- AMPLITUDE  SIMILARITY  SOLUTION 

The  equations  appropriate  to  large-amplitude  disturbances  admit  a  rich 
family  of  similarity  solutions,  which  illustrate  a  variety  of  behaviors. 
While  such  solutions  demand  particular  initial  conditions,  nonetheless  these 
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solutions  can  be  interpreted  as  good  local  approximations  for  situations 
arising  from  arbitrary  initial  conditions.   Two  cases  are  of  particular 
interest,  in  light  of  section  5:  one  describes  how  a  smooth  finite  peak 
evolves  to  an  infinite  singularity,  the  other  describes  how  that  first 
singularity  changes  shape  as  the  plume  evolves. 

Where  6  >>  de  (in  dimensionless  terms),  we  can  effectively  set  de  =  0  to 
get  the  large-amplitude  equations,  valid  locally.   In  fact,  the  same  equations 
would  apply  everywhere  if  somewhere  the  layer  has  zero  thickness;  for  the  case 
of  constant  end  forces,  minima  in  the  layer  (if  initially  less  than  de)  reach 
zero  thickness  in  finite  time.   Then: 

6t  +  (u6)x  =  0  (C.la) 

5  +  ux  =  0  (C.lb) 

We  assume  a  similarity  solution  of  the  following  form: 

x  =  a(t)C,  6  =  A(t)fU),  u  =  a(t)A(t)gU)  (C.2) 

Then  the  system  (C.l)  becomes: 

(A'/A2)f  -  (a'/(aA))£f  +  (fg)'  =  0  (C.3a) 

f  +  g'  =  0  (C.3b) 

Similarity  then  requires  that  A'/A2  be  a  constant.   If  that  constant  is 
not  zero  we  can  normalize  it  by  choosing  the  scale  of  A,  and  by  a  suitable 
choice  of  time  origin: 

A(t)  =  1/t  (C.Aa) 

Then,  because  6  must  be  non-negative,  solutions  where  f  >  0  will  apply  for 
t  >   0,  when  A  is  decreasing,  and  where  f  «  0  will  apply  for  t  <  0,  when  A  is 
increasing  catastrophically  to  t  =  0. 

Also,  a'/(aA)  must  be  constant,  say  X.   (As  X  is  insensitive  to  the  scale 
of  a(t),  it  cannot  be  normalized.)   Then  let: 
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a(t)  =  | t | A  (C.4b) 

For  t  >  0,  positive  X  gives  a  profile  that  spreads  out  in  the  x  direction, 

while  for  negative  X  the  profile  contracts;  the  converse  is  true  for  t  <  0. 
Combining  (C.3a,b): 

(gg')'  "  KZ"   -  g*  «  0  (C.5) 

which  integrates  to: 

(g-xe)g'  +  (X-l)g  -  C  (C.6) 

where  C  is  an  arbitrary  constant.   By  a  translation,  corresponding  to  moving 

the  origin  to  a  different  part  of  the  layer  moving  at  a  different  speed,  the 

constant  can  be  removed: 

?M-  C/(X2-X),  q(c)  =  g(0   ~   C/(X-1)  (C.7a) 

(q-Xc)q'  +  (X-l)q  =  0  (C.7b) 

as  long  as  \   #  0,1.   Then  the  equation  becomes  separable  for: 

Q(?)  e  cq(c) 

(Q-x)cQ'  +  (Q-DO  =  0  (C.8) 

Integrating  by  partial  fractions  and  substituting  back  gives: 

Z,   =  q  +  sgn(?-q)D|q|k  (C.9a) 

f  =  -1/[1  +  sgn(q(c-q))kD|q|k-l]  (C.9b) 

where  k  =  X/(X-1)  and  D  >  0  is  the  second  arbitrary  constant  of  integration. 
This  is  the  general  solution;  though  the  differential  equation  is  nonlinear,  a 
phase-plane  analysis  verifies  that  this  gives  all  solutions  (for  A(t)  not 
constant,  X  *   0,1)  except  the  trivial  solutions  f  =  0  and  f  =  -1. 

The  solution  describes  a  variety  of  behaviors,  depending  primarily  on  X 
and  on  which  branch  of  the  solution  is  chosen.  The  constants  C  and  D  affect 
the  x  origin  and  scale,  respectively. 
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For  example,  consider  \   >  1  (so  k  >  1),  C  =  0,  D  =  1,  and  the  particular 
branch: 

%   =  g  +  sgn(g)|g|k  (C.lOa) 

f  =  -1/[1  +  k|g|k-l]  (C.lOb) 

Since  f  <  0,  this  solution  only  applies  for  t  <  0;  it  describes  the  growth  of 
a  peak  up  to  the  singular  time  when  6  at  the  peak  becomes  infinite.   The 
asymptotics  in  £  reveal  the  behavior: 

as  |5|  +  0:  6  +  ! t  j  ~ 1 [ 1  -  k| 1 1 ~k | x |k"l ] ,  u  +  -x/ 1 1 1  (C.lOc) 

as  |^|  +  »:  6  -►  |x|"l/X/k,  u  +  -sgn(x)  | x |  l/k  (C.lOd) 

The  first  (C.lOc)  applies  for  x  near  the  origin,  as  long  as  t  is  non-zero. 
The  peak  is  of  fairly  general  shape,  but  for  the  peak  to  be  analytic  in  x, 
X  must  be  3/2  (k  =  3).   The  second  (C.lOd)  applies  far  from  the  origin  early 
on  (t  <<  0),  but  applies  ever  nearer  until  at  the  singular  time,  it  applies 
for  all  x  (#  0).   The  asymptotic  shape  (C.lOd)  is  independent  of  time;  as  the 
peak  grows  it  fills  in  the  integrable  negative  power  of  x  shape.   For  a  smooth 
peak,  X  =  3/2  and  at  the  singular  time  6  a    | x | —  2 / 3  (for  any  D) . 

As  another  example,  X  >  1,  C  =  sgn(^),  D  =  (X-l)k/X,  and  the  branch: 

£  =  g  -  sgn(g)X-l[((X-l)|g|  +  l)k  -  1]  (C.lla) 

f  =  l/[((X-l)|g|  +  l)k"l  -  1]  (C.llb) 

which  describes  a  plume,  symmetric  in  x,  for  t  >   0.   Asymptotically: 

as  |c|  +  0:  6  +  l/[  tl  -A/2/(  2  |x  |  )  ]  ,  u  +  -sgn(x)/(  2  |x  |  )/0~x^    (C.llc) 
as  |s|  +  oo:  6  +  (x|x|)"l/A,  u  +  -sgn(x)  (  X  |  x  |  )  1  /k/(X-l )       (C.lld) 
At  the  singular  time  (t  =  0),  the  profile  (C.lld),  which  is  independent  of 
time,  applies  everywhere.   In  other  words,  this  example  shows  what  happens  to 
an  initial  profile  proportional  to  an  integrable  negative  power  of  x;  by  a 
different  choice  of  C  and  D,  the  profile  (C.lld)  could  be  made  to  match 
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(C.lOd)  exactly.   For  t  >  0,  (C.llc)  applies  near  the  origin;  regardless  of 
the  starting  x  dependence,  the  shape  around  the  singularity  (plume)  becomes 
proportional  to  l//|x|.   If  X  <  2  then  the  square-root  singularity  and  the 
corresponding  strength  of  the  sink  at  the  origin  decay  with  time,  as  the  layer 
drains  away.   (This  would  be  true  for  an  initially  smooth  profile  like  (C.lOc) 
with  X  =  3/2.)   Conversely,  for  X  >  2,  the  singularity  grows,  as  fluid  comes 
in  from  the  sides  faster  than  it  can  be  disposed  of.   The  special  case  X  =  2 
gives  a  steady  plume: 

6  =  l//(2|x|),  u  =  -sgn(x)/(2|x|)  (C.lle) 

Also,  there  are  related  solutions  for  X=0,  0  <  X  <  1,  and  X  =  1,  which  have 
the  same  asymptotic  shapes  (C.llc,d),  except  that,  for  |  £  |  ■*■  °°,  the  forms  for 
u  are  different,  and  for  X  =  0,  6  ■*■  e~lxl/t  as  |  £  |  +  °°. 

Briefly,  the  other  behaviors  governed  by  the  similarity  solution  are  as 
follows.   For  t  <  0,  i.e.,  profiles  growing  to  the  singular  time,  there  are 
four:  a  minimum  thickness  of  zero  like  |x|^,  k  >  0,  locally  time-independent, 
that  far  away  levels  off  to  approach  l/|t|;  a  profile  that  approaches  zero  as 
x  -►  -oo  and  l/|t|  as  x  -)-  +»;  a  symmetric  finite  minimum  flanked  by  plumes 
(which  could  be  extended  periodically);  and  a  plume  whose  sides  level  off  to 
approach  l/|t|  (rather  than  0).   For  t  >  0,  where  the  profile  starts  at  the 
singular  time  and  diminishes  thereafter,  there  is  only  one  other  case:  a  zero 
minimum  flanked  by  plumes  (which  could  extend  periodically).   The  special  case 
of  A(t)  =  1 ,  so  a(t)  =  e^t,  gives  either  an  isolated  plume,  or  a  zero  minimum 
where  |  6X  |  ->■  °°  surrounded  by  plumes  (possibly  periodic).   All  cases  with 
plumes  show  the  l//|x|  shape.   And  solutions  with  the  same  A(t)  and  X  but 
different  C  and  D  can  be  "cut  and  pasted"  together  to  give  other  similarity 
profiles;  if  g  and  g"   are  continuous  this  will  produce  reasonable  results. 
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Figure  1.   Raylelgh-Taylor  Problem.   A  layer  of  fluid  1  over  a  layer 
of  a  different  fluid  2,  between  shear-free  boundaries,  where  Po  <  Pi  • 
Interface  position  given  by  z  =  6(x,t).   Equilibrium  (unstable;  depths 
of  upper  and  lower  layers  are  dj  and  d2«   Growth  rate  of  instability 
depends  on  viscosity  ratio  a  =  uj/u2>  thickness  ratio  8  =  dj/d2,  and 
wavenumber  k  (dimensionless  e  =  kdj). 


Figure  2.   Viscous  Layer  Profiles  6(x)  for  various  times  t,  from 
initial  sinusoidal  perturbation  of  amplitude  10~4,  for  constant  end- 
force  conditions  (de  =1).   Each  profile  represents  one  wavelength; 
the  overall  strain  increases  with  time. 
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Figure  3.   Viscous  Layer  Profiles  <5(x)  for  various  times  t,  from 

initial  sinusoidal  perturbation  of  amplitude  10"^,  for  constant 
wavelength  conditions  (de(t)  decreases  with  time).   Similar  to  Fig.  2, 

except  troughs  not  as  deep,  and  the  peak  reaches  °°  slightly  sooner. 
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Figure  4.   Peak  Disturbance  Growth:  «max(t)-l  (logarithmic)  vs.  time, 
from  initial  sinusoidal  perturbation  of  amplitude  10-4. 

(solid):  for  constant  end  force 

(dashed):  for  constant  wavelength 

(dotted):  exponential  (linearized)  growth,  for  comparison 
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Figure  5.   When  small-slope  approximation  fails,  it  does  so  only  in  a 
region  of  width  0(e) . 
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Figure  6a.   Linearized  Growth  Rate:  o(k)  (logarithmic  scales) 

(a)  Infinite  depth  (&+">),    for  various  viscosity  ratios  a.   This 
shows  that  when  the  film  is  much  more  viscous  than  the  fluid  below 
(a<<l),  the  maximum  growth  rate  applies  over  a  broad  range  of 
wavelengths  (to  lowest  order). 
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Figure  6b.   Linearized  Growth  Rate:  a(k)  (logarithmic  scales) 

(b)   Finite  depth  (various  6),  for  a  =  10~5.   This  shows  how 
finite-depth  effects  give  a  large-wavelength  cutoff  to  the  maximum 
growth  regime  for  o«i. 
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Figure  7.   Rayleigh-Taylor  Growth  Regimes.   The  various  (dimension- 
less)  asymptotic  growth  rates  o  (from  Eq.  k._2>c)   are  shown  in  the 
regions  of  the  parameter  plane  (wavenumber  k  and  viscosity  ratio  a, 
logarithmic  scales)  where  they  apply,  assuming  the  depth  ratio  6>>1. 
(Divisions  between  regions,  shown  as  lines,  are  actually  broad  areas 
across  which  the  growth  rate  varies  continuously.)   These  growth 
regimes  apply  while  the  interface  slope  remains  small,  which  for  lonj 
waves  (k<<l)  includes  large  amplitudes. 
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Figure  8a.   Power-Law  Fluid  Layer  Profiles:  6(x)  for  various  times  t, 
from  an  initial  sinusoidal  perturbation  of  amplitude  0.02,  assuming 
constant  end-force  conditions, 
(a)   for  power  m  «=  3 

(dotted):  initial,  t  =  0 

(dashed):  t  =  129 

(solid):  t  =  153 
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Figure  8b.   Power-Law  Fluid  Layer  Profiles:  6(x)  for  various  times  t, 
from  an  initial  sinusoidal  perturbation  of  amplitude  0.02,  assuming 
constant  end-force  conditions, 
(b)   for  power  m  *  9 

(dotted):  initial,  t  =  0 

(solid):  t  =  8.8*109 
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Figure  9a.   Power-Law  Fluid  Disturbance  Growth:  6(t)  (linear  scales) 
following  a  fluid  cross  section,  from  an  initial  perturbation 
amplitude  of  0.02,  assuming  constant  end  forces.   Growth  is  shown  for 
both  positive  and  negative  perturbations, 
(a)   for  power  m  =  3 
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Figure  9b.   Power-Law  Fluid  Disturbance  Growth:  <5(t)  (linear  scales) 
following  a  fluid  cross  section,  from  an  initial  perturbation 
amplitude  of  0.02,  assuming  constant  end  forces.   Growth  is  shown  for 
both  positive  and  negative  perturbations, 
(b)   for  power  m  =  9 
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Figure  10.   Force  balance  on  a  differential  column  of  fluid  (z 
upward) . 
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